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ATTITUDE/ATTITUDE-RATE ESTIMATION 
FROM GPS DIFFERENTIAL PHASE MEASUREMENTS 
USING INTEGRATED-RATE PARAMETERS 

Yaakov Oshman* and F. Landis Markley t 


A sequential filtering algorithm is presented for attitude and attitude-rate estima- 
tion from Global Positioning System (GPS) differential carrier phase measurements. 
A third-order, minimal-parameter method for solving the attitude matrix kinematic 
equation is used to parameterize the filter’s state, which renders the resulting estima- 
tor computationally efficient. Borrowing from tracking theory concepts, the angular 
acceleration is modeled as an exponentially autocorrelated stochastic process, thus 
avoiding the use of the uncertain spacecraft dynamic model. The new formulation 
facilitates the use of aiding vector observations in a unified filtering algorithm, which 
can enhance the method’s robustness and accuracy. Numerical examples are used to 
demonstrate the performance of the method. 


INTRODUCTION 

Attitude determination methods using Global Positioning System (GPS) signals have been inten- 
sively investigated in recent years. In general, these methods can be classified into two main classes. 
Point estimation algorithms (also called “deterministic” algorithms), in which the GPS measure- 
ments at each time point are utilized to obtain an attitude solution independently of the solutions at 
other time points, were introduced, among others, in Refs. 1, 2 and 3. Stochastic filtering algorithms, 
which process the measurements sequentially and retain the information content of past measure- 
ments, can produce better attitude solutions by more effectively filtering the noisy measurements. 
Such algorithms were recently introduced in Refs. 4 and 5, both of which utilized extended Kalman 
filtering to sequentially estimate the attitude from GPS carrier phase difference measurements. Both 
attitude and attitude- rate were estimated, and the filters used the nonlinear Euler equations of mo- 
tion for attitude propagation. While avoiding the traditional usage of the costly and unreliable gyro 
package, this approach rendered the resulting filters computationally burdensome and sensitive to 
inevitable modeling errors. 6 In Ref. 4 an attempt was made to robustify the dynamics- based filter 
by estimating the unknown disturbance torques, modeled as unknown constants. 

Although GPS-based attitude estimation methods should enjoy, in principle, the low price and 
low power consumption of state-of-the-art GPS receivers, and the general availability and robustness 
of the global positioning system, these methods are very sensitive to multipath effects and to the 
geometry of the antennae baseline configuration, and they inherently rely on precise knowledge 
of the antennae baselines in the spacecraft body frame. On the other hand, methods based on 
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vector observations have reached maturity and popularity in the last three decades. However, as is 
well known, they too suffer from disadvantages, that can be attributed to the particular attitude 
sensors on which they are based. Thus, while their readings are relatively noiseless, Sun sensors are 
very sensitive to Earth radiation effects, and are rendered completely useless during Eclipse. Star 
trackers can provide accuracy on the order of a few arc-sec, but are usually extremely expensive. 
Magnetometers always provide measurements of the Earth magnetic field in spacecraft flying in low 
Earth orbits, but they are sensitive to unmodeled residual magnetic fields in the spacecraft and to 
magnetic field model imperfections and variations. 

The method presented herein is a sequential estimator for both the spacecraft attitude matrix and 
attitude-rate, which mainly uses differential GPS carrier phase measurements, but can also process 
aiding vector observations (such as low accuracy coarse Sun sensor measurements, or magnetic field 
measurements). Conceptually similar to the principle of complementary filtering, the idea underlying 
this estimator is that, due to the different nature of these signals, the combination of both in a unifi ed 
data processing algorithm can benefit from the relative advantages of both sensor systems, while 
alleviating the disadvantages of both. 

The new estimator is based on a third-order minimal-parameter method for solving the attitude 
matrix evolution equation using integrated-rate parameters (IRP). 7 Similarly to Refe. 5 and 4, 
the new estimator is a sequential filtering algorithm and not a deterministic (point estimation) 
algorithm. However, the new algorithm differs from other works addressing the same problem in 
two main respects. First, the estimator’s propagation model does not utilize the nonlinear Euler 
equations. Instead, employing an approach borrowed from linear tracking theory, 8 the uncertain 
dynamic model of the spacecraft is abandoned, and the angular acceleration is modeled as a zero- 
mean stochastic process with exponential autocorrelation. Combined with the extremely simple 
evolution equation of the integrated-rate parameters, this results in a simple, linear propagation 
model. Second, in contrast with other methods relying m ainl y on the attitude quaternion, the 
algorithm presented herein directly estimates the attitude matrix, a natural, nonsingular attitude 
representation. Building upon the min i m al, third-order integrated-rate parametrization, the new 
estimator assigns just three state variables for the parametrization of the nine-parameter attitude 
matrix, which is at the heart of its computational efficiency. 

INTEGRATED-RATE PARAMETERS 

Consider the matrix differential equation 

V(t) = W(t)V(t), V(to) = Vo (1) 

where V(t) 6 R n - n , W{t) = -W T (t) for all t > t 0> V 0 V? = I and the overdot indicates the temporal 
derivative. De fining 


A{*,to)± fw{r)dr (2) 

Jto 

W Q (t)±W(t)-(t-to)W(t) (3) 

it can be shown that the following matrix-valued function is a third-order approximation of V'(f): 

V(t,t 0 ) = {/ + A(t,t 0 ) + + A3 (*j*o) +Lzlo[ A(t>to) fv 0 (t) - W 0 (t)A(t,t 0 )] }V 0 (4) 

Moreover, V is a third-order approximation of an orthogonal matrix, i.e., V(t,t 0 )V T (t,t 0 ) = I + 
0((f - t 0 ) 4 ) where 0(x) denotes a function of x that has the property that 0{x)/x is bounded as 
x -» 0. 

In the 3-D case, the off-diagonal entries of -4(<,t 0 ), termed integrated-rate parameters, have a 
simple geometric interpretations they are the angles resulting from a temporal-integration of the 
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three components of the angular velocity vector c o(t) = [u;i(i) W 2 (t) uz(t)] T , where t Ji is the 
angular velocity component along the i-axis of the initial coordinate system, and i = 1,2,3 for 
x,y, z, respectively. The orthogonal matrix differential equation (1) is rewritten, in this case, as 

D(t) = n(t)D(t), D(t 0 ) = D 0 (5) 

where D{t) is the attitude matrix, or the direction cosine matrix (DCM), f l(t) = — [u;(t)x], and 
[u;(t)x] is the usual cross product matrix corresponding to u;(t). In this case, the matrix A(t,to) 
takes the form 


>l(Mo) = -[*(*) *] 

where the parameter vector d(t) is defined as 

«(*)±[*l(«) 02 (t) 

and 


9i(t)= [ Ui{T)dT , 

Jt o 


*« 1.2,3 


( 6 ) 


( 7 ) 


( 8 ) 


Let the sampling period be denoted by T ^ t^+i - t*. Using the notation 9{k) = 0(th)> the 
parameter vector at time t * is 9(k) = [Oi(k) 02{k) ^3(^)] r and Eq. (8) implies 


f tk 

&i(k) = / Wi(r)dr, * = 1,2,3 
Jt 0 

From Eq. (9) we have 

r*k + i 

9(k -f 1) = 9(k) 4* / u>(r) dr 

Define A(k 4* 1, Ar) to be the discrete-time analog of A(t,to), he., 

i4(fc + 1, Jfc) 4 - [(*(* + l) - 9(k)) x] 

Also, let >F(A: + 1) = — [V>(& 4- l)x], where 

\l>(k 4- 1) — u;(A; + 1) — clj(k 4- 1)T 
Then, the corresponding discrete-time equivalent of Eq. (4) is 

D(k + 1) = j J + A(k + 1, jfc) + ^A*(k + 1, Jfc) + i A 3 (k + 1, k) 

+ |t [A(k + 1, J fc)tf (jfc + 1) - *(k + l)A(jfc + 1, jfc)] |d(A) 

which, using Eqs. (11) and (12), can be written as 

D(k + 1) = D[d[k + 1) - 0(jfc),w(jfc + 1 ),w(Jb + l),D{k)] 


( 9 ) 


( 10 ) 


( 11 ) 


( 12 ) 


(13) 


(14) 
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KINEMATIC MOTION MODEL 

To avoid using the uncertain spacecraft dynamic model, the spacecraft angular acceleration is mod- 
eled as a zero-mean stochastic process with exponential autocorrelation function. The acceleration 
dynamic model is, therefore, the following first-order Markov process, 

o?(t) = -Ad;(t) + i>(t) (15) 

For simplicity, a decoupled kinematic model is chosen for the three angular rate components, i.e., 
A = diagjrj 1 ,r 2 1 ,r 3 *}, where {ri} 8 =1 are the acceleration decorrelation times associated with 
the corresponding body axes. The driving noise is a zero-mean white process, with power spectral 
density (PSD) matrix 


Q(t) = 2AE 2 , E = diag{<7i, 0 * 2 , cr 3 } (16) 

The noise variances in Eq. (16) were chosen according to the Singer angular acceleration probabilis- 
tic model, 8 in which the angular acceleration components, can be 1) equal to cJMi with 

probability p Mi1 2) equal to -Co Mi with probability p Mi , 3) equal to zero with probability po», or 
4) uniformly distributed over the interval [— wjv/i,cjjvf *] with the rem aining probability mass. Using 
this model, it follows that 


= — ^p-O- + 4PAfi - POi) (17) 

The parameters Cjm pM i and po* are considered as filter tuning parameters. As customarily done, 
they are selected by experience with real and simulated data, so as to optimally adapt the filter to 
the characteristics of the problem at hand. 

Now let the system's state vector be defined as x(t) = [0 r (t) v T (t) u T {t )) T , then the state 
equation is 



0 1 O' 


r o i 

±{t) = Fx(t) + v(t) = 

0 0/ 
0 0 -A 

x(t ) + 

* 

3 ° 


(18) 


with obvious definitions of F and v(t). Corresponding to the sampling interval T, the discrete-time 
state equation is 


x(k + l) = 9(T)x(k)+v(k) 


(19) 


where the transition matrix is 


$(r) = e FT = 


/ 

0 

0 


TI A -2 (e~ AT — I + TA) 
I A -1 (J - e -AT ) 

0 e‘ AT 


and v(k) is a zero-mean, white noise sequence, with covariance matrix 


Q(k) A E{v(k)v T (k)} 



e F < r -‘>diag{0,0 ,Q(t)}e FT ^dt 


( 20 ) 


( 21 ) 


MEASUREMENT PROCESSING 
GPS Differential Phase Measurements 

Consider the basic GPS antenna array, depicted in Fig. 1. The array consists of the master ante nna, 
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Figure 1. GPS Phase Difference Measurement Geometry 


A m , and the slave antenna, Aj. These antennas are located on the satellite’s surface, such that the 
baseline vector between them, resolved in a body-fixed coordinate system, is bj. It is assumed that 
the entire system consists of m*, antennas, in addition to the master antenna, so that there exist m*, 
independent baselines. It is also assumed that at time GPS satellites are in view. 

Consider the ith satellite, and denote the sightline (unit) direction vector to that satellite, resolved 
in an inertial coordinate system, by s*. Let D(k 4 1) be the attitude matrix transforming vectors in 
the inertial coordinate system to their body-fixed system representations at time Let Nij(k+1) 
and A&j(A; + l) denote the integer and fractional parts, respectively, of the phase difference between 
the two carrier signals, corresponding to the ith satellite, as acquired by the antennas A m and Aj. 
Denoting by A the GPS carrier wavelength, the true (noiseless) signals satisfy 

[A 4>ij{k + 1) + Nij(k + 1)]A = bjD(k + l) Si (22) 

The standard GPS carrier wavelength is 19.03 cm. In this work, it is assumed that the integer part 
of the phase difference between the two receivers is known from a previous solution. 1,9 

In practice, the phase measurements will be contaminated by noise, the primary source of which 
is due to the multipath effect. 1 Denoting the noise corresponding to the baseline bj and the sightline 
Si by fiij(k -hi), the read measurement equation is 

[A <j>ij(k 4 1 ) + Nij(k + 1)] A = bjD{k 4 l)si 4 hij(k 4 1) (23) 

where it is assumed that hij{k 4 1) ~ N(0 , ^(k -h 1)). Typically it can be assumed that the noise 
standard deviation is on the order of 5 mm. 1 From Eq. (23) we obtain the normalized measurement 
equation 


A + 1) + Nij(k 4 1) = bjD(k 4 1 )si 4 riij(k 4 1) (24) 

where we have defined bj = bj /A and n t j(k -h 1) = fiij{k -h 1)/A. The normalized measurement noise 
satisfies riij(k + 1) ~ N(0, cr?j(k + 1)), where a x j(k + 1) = &ij(k + 1)/A. 

GPS Measurement Linearization 

At tk+i the minimum meam square error (MMSE) predicted vector is x(fc+l|A;), and its corresponding 
prediction error covariance matrix is P(k 4- l|/c) = E{x(k 4 l|A:)f r (A: 4 1| fc)}, where the estimation 
error is x{j\k) = x(j) - x(j\k). Using Eq. (14), Eq. (24) is rewritten as 

Nij(k 4 1) 4 A <(>ij(k 4 1) = bjD[9(k 4 1) - 6{k),u{k 4 1 ),u>(fc 4 1), D(/c)]si 4 riij{k 4 1) (25) 
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Next, we linearize the nonlinear measurement equation (25) about the most recent estimate at tk+i, 
i.e., 


x(k + 1) = x(k + 1|A:) + Sx(k + 1) = 


~9(k 4- 1|A:)" 


'89(k + 1)' 

u(k 4- l|fc) 


6cj(k H- 1) 

jj(k + l|fc)_ 


Suj(k 4- 1) 


(26) 


where 89(k + 1), 8u>(k + 1) and 5u(k + 1) are the perturbations of the state components about the 
nominal (i.e., predicted) state. Let D*(k\k) denote the a posteriori, orthogonalized estimate of the 
attitude matrix at time f*, to be discussed in the next section. Using now the most recent estimates 
for D(k ) and x(k), namely D*(k\k) and x(/c|&), respectively, in Eq. (25), it follows that 

A< M fc + 1) + N v (k + 1) = bj D[9{k + l|ft) + 89(k + 1) - ^(Ar|fc),d>(A: + l\k) + 5u(k + 1), 

u(k + 1|A) + 6u(k + 1), D* (*!*)]«< + mj(k + 1) (27) 

As discussed in the sequel, the a posteriori IRP estimate is zeroed after each measurement update 
(due to full reset control of the IRP state). We will, therefore, use the reset value of the IRP estimate, 

9 (fc|&) = 0, in Eq. (27). Now expand D about the nominal state using a first-order Taylor series 
expansion, i.e., 


D[9{k + 1|A) + S9(k + l),u(k + l|fc) + 5w{ k + l),<L(k + 1|*) + 8<j(k + 1), £>*(Jfc|Jt)] 

59i(k + 1) 


= D(k + 1| k) + V^ g ~ P [ g ( /: + 1 )' (2, ( fc + 1 l fe ).^(fe+l|fe),-P*(fc|A:)] j 

99 i 




dD\9(k + l|fc),o>(A: + l),d>(fc + ![/;), D»(fc|fc)] 


t=l 

3 


dcji 


w(*+l|*) 


+ 


dD[§(k 4* l|fc) 7 C){k + l|t),o;(Ar 4- l),J9*(A:|fc)] 

^ dUi 


i= 1 


w(*+l|*) 


0(fc+l|fc) 
6uJi(k -f 1) 

Sui(k -I- 1) 


(28) 


where (»)|^ denotes ‘evaluated at C’ and £>(fc + l\k) A Z)[»(* + l|Jfc),w(A: + l|A;),J;(fc + l|A),D*(Jfc|A;)] 
Differentiating Eq. (13), the sensitivity matrices appearing in Eq. (28) are computed as 
d 

OfD[9{k + l),w(A: + l|Ar),£^(A: + l|fc), £)*(fc|*r)] = <?i[0(A: + 1),^(A; + l|Jfc)]D*(jfc|A:) (29a) 

d 

^D[0(fc + l|A:),a;(/: + l),^ + l|fc),£)*(*|fc)] = ^TFi [9(k + l|Jk)]D*(A;| k) (29b) 

^D[9{k + l|fc),w(fc + 1|A:), w(fc + 1), D*(fc|fc)] = -br 2 F t [9{k + l|fc)].D*(Jfc|/fc) (29c) 

for t = 1,2, 3, where tp(k + 1|A:) = u;(k + l|fc) - Tu)(k + l|Jk), and 

Gi{9, rp) = -(9e[ + erf 7 ) - 9J - (1 - i ||^|| 2 ) [e*x] + gT(^ef - e i V’ T ) + ^[flx] (30a) 

F i (9) = e i 9 T -9e 7 (30b) 

where e, is the unit vector on the zth axis, i = 1,2,3. 

Using Eqs. (28), (29) and (30) in Eq. (27) yields 

A<t>ij(k + 1) + Nij(k + 1) - bj£)(k + l|fc)sj = hj^k + l)8x(k + 1) + ny(Jfc + 1) (31) 
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where the observation vector hij(k 4 1) € R 9 is defined as 

M* + l)= [*•«(* + 1) ^«(* + l) hjjik + l)] 7, (32) 

and the elements of the vectors he i j(k 4- 1) <E K 3 , h UX j{k -f 1) € R 3 and h i j i j(k -f 1) € R 3 are 


h Oij p (k + 1) = bjG p [9(k 4- 1| k),i>(k + l|fc)]£»*(fc|A:)s I , 

p = 1,2,3 

(33a) 

h uijp {k + 1) = ±TbjF p [6(k + 1| k)]D*(k\k) Si , 

P = 1,2,3 

(33b) 

hujijpik 4 1) = — Th ux jp(k 4 1), 

P = 1,2,3 

(33c) 


Define now the effective GPS measurement to be 

ytj(k + 1) 4 Afrjik + 1) + N tj - bjD(k + l|fc) 5< (34) 

Then, using this definition in Eq. (31) yields the following scalar measurement equation: 

yf.(k + 1) = hi/ik + 1 )Sx(k + 1) + 7 Uj(k + 1) (35) 

For the m*> baselines and m 3 sightlines, there exist m 3 x scalar measurements like Eq. (35). 
We next aggregate all of these equations into a single vector equation, such that the measurement 
associated with the baseline bj and sightline $i corresponds to the pth component of the vector 
measurement equation, where p = (j — l)m s 4 i. This yields 

y*(k + 1) = H*(k 4 1 )5x(k 4 1) + n*{k 4 1) (36) 

where thepth row of the matrix H^(k-\-l) is hij T (k + 1), n^(A;4l) ~ N(0, R?(k+1)), and i?^(A;4l) 
is a diagonal matrix whose diagonal elements are R£p(k 4 1) = a\y 


Vector Observation Aiding 

If the sole source of attitude information is the GPS carrier phase signals, then Eq. (36) should serve 
as the basis for the development of the measurement update algorithm (in the next section). In the 
case that vector observations are available, this information structure needs to be augmented. 

Assume that a new pair of corresponding noisy vector measurements is acquired at tk+i* This 
pair consists of the unit vectors u(k 4 1) and v(k + 1), which represent the values of the same vector 
r(A;4l), as modeled in the reference coordinate system and measured in the body coordinate system, 
respectively. The direction-cosine matrix D(A:4 1) transforms the true vector representation uq into 
its corresponding true representation vo according to 


vo (k 4 1) = D(k 4 l)uo(fc 4- 1) 


(37) 


Assuming no constraint on the measurement noise direction, the body-frame measured unit vector, 
v(k 4- 1), is related to the true vector according to 


_ vq ( k 4- 1) 4* n f v (k -I- 1) 

IM* + 1) + n v(k + 1)11 


(38) 


where the white sensor measurement noise is n' v (k 4- 1) ~ N(0, R' v (k -I- 1)). Since both vo(k 4- 1) and 
v(k 4 1) are unit vectors, it follows from Eq. (38) that 


v(k 4 1) = vo (k 4 1) 4- n v (k 4 1) 


(39) 
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where n v (k 4 1) — 7y 0 (k 4- l)n' v (k 4 1) and 7„ 0 (k 4- 1) — I — vo(k 4* IJvq’ (A: 4 1). To a good approx- 
imation, the effective measurement noise is a zero mean, white Gaussian sequence with covariance 

Rv(k 4 1) = ?i(k 4 1 )R' v (k 4 l)?t 0 (k 4 1) (40) 

To account for non-ideal effects (e.g., star catalog errors), it is assumed that the modeled reference 
vector is related to the true vector according to 

u(k 4 1) = uo(A: 4 1)4 ti u ( k 4 1) (41) 

where J_ uo is a zero mean, white Gaussian noise, that is uncorrelated with n v and has a known 
covariance matrix Ru(k). 

Vector Measurement Linearization 

Using Eqs. (11), (12) and (13), Eq. (37) can be rewritten as 

v 0 (k + 1) = D[9(k + 1) - 9(k),w(k + 1 ),u(k + 1), D(k)]uo(k + 1) (42) 

Linearizing about the predicted estimates and using Eqs. (26), (39) and (41), it follows that 

v(k + 1) - ^{k + 1) = D[9{k + l|fc) + 69(k + 1), u(k + l|Jb) + 5w{k + 1), 

u(k + l|fc) + 5u{k + 1), £>*(A|k)] [u(k + 1) - n u (k + 1)] (43) 

where the reset value of the IRP estimate, 0 c (k|k) = 0, has been used. Expanding D about the 
nominal state using the first-order Taylor series (28) yields 

3 

v(fc + 1) - D(k + l\k)u(k + 1) = ^ [°i [*(* + 1|*)> *(* + l|*)]^i(Jk + 1) 

t=l 

+ ^Fi[0(A + 1|*)]&*(* + 1) - 1 T 2 Fi[9(k + l|fe)]^i(fc + l)]£>*(Jfc|*)u(A: + 1) 

- D(k + 1|*K(* + 1) + n v (k + 1) = H u {k + l)Sx(k + 1) - D{k + l|*)n„(Jfe + 1) + n»(lfc + 1) 


(44) 

where the observation matrix H u {k + 1) is written in block matrix form as 

H v (k + 1) = [fli (k + 1) H 2 (k + 1) H 3 (k + 1)] € R 3 - 9 (45) 

and the columns of the submatrices Hi(k + 1) € R 3,3 , i = 1, 2, 3 are 

Hij(k + 1) — Gj [9(k + l|k),^(A: + l|fc)].D*(A:|A:)u(A: + 1) (46a) 

H 2j (k + 1) = ±TFj [9(k + l|*)]£-(Jb|Jfc)t*(Jfe + l) ( 46b ) 

H 3 j(k + 1) = -TH 2j (k + 1) (46 C ) 

^ 3 L 2, 3. Define now the effective measurement and measurement noise to be, respectively, 

y u {k + 1) 4 v [k + 1) - D(k + l|Jfc)u(jfc + 1) (47) 

n v {k + 1)4 ^(k + 1) - D(k + l(Jb)n„(A: + 1) (48) 
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Then, using these definitions in Eq. (44) yields the following measurement equation: 

y u (k + 1) = H u {k + 1 )Sx(k + 1) + n u {k 4 1) (49) 

where n v (k -h 1) ~ N(0, R u (k 4- 1)) is the white measurement noise, and 

R u (k + 1) = Rvik + 1) + D(k + ll^Rvik + l)£> r (ifc + 1| k) (50) 

Measurement Update 

To process the measurements, define now 



where n ~ 3sf(0 y R) and R = diagl^, ii 17 }. Since Sx(k 4* 1) = x(k 4- 1) — x(k 4- 1| k) = x(k + l|fc) 
and x(k 4- 1| k) is an unbiased, MMSE predictor, we have E{6x(k + 1)} = E{x(k 4- 1 1 A:) } = 0 and 
cov{6x(k 4* 1)} = cov{x(A 4- l|fc)} = P(k 4- 1| A), thus 5x(k 4- 1) ~ X(0,P(A 4- 1| A)). Using the 
linearized measurement equation and the statistical properties of the measurement and prediction 
errors, the MMSE estimator of Sx(k 4- 1) is 

Sx(k 4- 1|A 4- 1) = K{k 4- l)j/(A 4- 1) (52) 

where K(k 4- 1), the estimator gain matrix, is computed as 

K(k + 1) = P(k + l\k)H T (k + 1) + 1)P(& + l|Jfc)fT T (jfc + 1) + R(k + 1)] _1 (53) 

Also, Jx(A+l|A4-l) = x(A4-l|A4-l)-x(A4-l|A) which, used in Eq. (52), yields the state measurement 
update equation 


x(A 4- 1 1 A: 4- 1) = x(A 4- 1|A) 4- K(k 4- 1 )y{k 4- 1) (54) 

Subtracting x(A 4- 1) from both sides of the last equation yields 

x(A 4- 1| k + 1) = [I - K(k 4- 1 )H(k 4- 1 )]x(A + l|Ar) - K(k + l)n(A 4- 1) (55) 

from which the resulting covariance update equation is 

P{k + l|Ar + 1) = [/ - K(k + l)H(k + 1)] P(k + l|Jk) [/ - K(k + \)H{k + 1)] T 

+ K(k + l)R(k + l)K T (k + 1) (56) 

where the filtering error covariance is P(k + 1 [A: + 1) = E{x(k + l|fc + l)x T (k + l|fc 4- 1)}. 

To compute the measurement-updated attitude matrix at time t^+i, we use the most recent 
estimate x(A + 1|A4- 1) and the estimated attitude matrix corresponding to time tk in Eq. (13). This 
yields 


D(k + l|fc + 1) = {I + A(k + l,k) + lA 2 {k + l,k) + ^A 3 (k + l,k) 


4 


4- 


ir[i(fc + 1, k)*(k + 1| k + 1) - 4f(k + l|Jfc + l)i(A: 4- 1, Jfe)] jl>*(A:|Jfc) (57) 
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where the a posteriori estimates of A(k + 1, k) and ^(k 4* 1) are defined, respectively, as 

i4(* + l,fc)A-($(* + i|* + i)x] t $(Jfe + l|* + l)±-[tf(* + l|* + l)x] (58) 

where V>(/c + l|/: + l) = w(A;+l|/i: + l) — Tw(fc+l|/: + l), and D*(k\k) is the a posteriori, orthogonalized. 
estimate of the attitude matrix at time tk, to be discussed in the next section. 

Finally, since the a posteriori attitude matrix, D{k + l|ifc + 1), is computed based on the a 
posteriori estimate, 9{k + 1 1 Ar + 1), this implies a full reset control of the parameter vector, i.e., 
9 c (k + 1) = 9(k + 1) — 6(k + 1|A: 4- 1), where 0 c (k + 1) is the reset state vector at tk+ 1 , and a 
corresponding reset of the state estimate, 6 c (k + 1 |Ar + 1) =0, which is then used in the ensuing 
time propagation step. Since the reset control is applied to both the state vector and its estimate, 
no changes are necessary in the estimation error covariance matrix. 

ATTITUDE MATRIX ORTHOGONALIZATION 

To improve the algorithm’s accuracy and enhance its stability, am additional orthogonalization pro- 
cedure is introduced into the estimator, following the measurement update stage. In this procedure, 
the orthogonal matrix closest to the filtered attitude matrix is computed. 

Given the filtered attitude matrix D(k + l|fc 4- 1), the matrix orthogonalization problem is to 
find the matrix 

D*(k + l|fc + 1) 4 arg^ncdn 3 ||l>(A; 4- l|fc + 1) - £>|| , subject to D T D = I (59) 

Being a special case of the orthogonal Procrustes problem, the matrix orthogonalization problem 
can be easily solved using the singular value decomposition (SVD). In cases where the excessive 
computational burden associated with the SVD might render its use prohibitive, e.g., in real-time 
attitude determination and control, the following approximate orthogonalization method, based on 
the iterative method introduced in Ref. 10, cam be utilized: 

D*(k + l\k + l) = N(k + 1 )D(k + l|Jb + 1) (60) 


where 


N(k + 1) A |/ - lb[k + 1|* + 1 )D T (k + l|fc + 1) (61) 

Remark 1. Using an approach similar to that used in Ref. 11, it can be shown that, to first-order 
accuracy, the orthogonalization procedure does not affect the statistical properties of the estimator 
and, therefore, does not necessitate any adjustments in the algorithm. 

PREDICTION 

In the prediction step at , the reset a posteriori estimate at time tk, x c (k\k) (computed with the 
reset IRP estimate) and its corresponding error covariance matrix, P(k\k), are propagated to tim e 
lk + 1 ■ 

Using Eq. (19), we have 


x(k + 1|A?) = $(T)x c (Jfc|jfc) 

Using this result with Eq. (19) yields the covariance propagation equation 
P(k + 1| k) = $(T)P(k\k)$ T (T) + T(T)Q f (k)r T (T) 


(62) 


(63) 
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To propagate the attitude matrix to tk+i we use the most recent IRP, attitude-rate and angular 
acceleration estimates, and the orthogonalized DCM estimate corresponding to £*, in Eq. (13). This 
yields 


D(k + l\k) = ll + A(k + l,k) + lA 2 (k + l,k) + lA 3 (k + l,k) 

( 2 6 

+ ^T^A(k + 1, k)V(k + l|fc) - 4>(k + l|Jfc)J(ifc + 1, *)] |ir(Jfc|Jk) 

where the a priori estimates of A(k + 1, k) and ^(/c + 1) are defined, respectively, as 
A(k + l,k) 4 — [^(Jfc + l|fc)x], ^ -[^(Jfc+l|ifc)x] 


(64) 

(65) 


NUMERICAL STUDY 
Example I 

In this example, three non-orthogonal baselines were used: bi = [l.O, 1.0, 0.0] T , 62 = [0.0, 1.0, 0.0] r , 
63 = [0.0, 0.0, 1.0] T . Two fixed sightlines were observed at all times, si = -^[l.O, 1.0, 1.0] T and 

S2 = ^-[0.0, 1*0, 1.0] T The non-normalized GPS signal noise standard deviation was 5.0 mm. 
When vector measurements were used, the noise equivalent angle of the inertially- referenced obser- 
vations was set to 5.0 arc-s, while the body- referenced vector measurements were simulated to be 
acquired by a low accuracy attitude sensor with a noise equivalent angle of 0.1 deg. These mea- 
surements corresponded to a randomly selected vector, which was kept constant throughout the 
run. 

The angular rates of the satellite satisfied u ;*(t) = Aisin^t + <fo), where A{ = 0.02,0.05,0.03 
deg/s, 4>i = 7t/4, 7r/2, 37r/4 rad, and 7i = 85,45,65 s for i = 1,2,3, respectively. The initial angular 
rate estimates were all set to zero. The true initial attitude corresponded to Euler angles of 30 deg, 
20 deg and 10 deg in roll, pitch and yaw, respectively, while the filter’s initial state corresponded to 
Euler angles of 25 deg, 15 deg and 5 deg, respectively. The filter was run at a rate of 20 Hz, and 
the measurement processing rate was 10 Hz. The Singer angular acceleration model was used with 
parameters set to r = 10 s, Cjm = 10 -4 rad/s 2 , pm = po = .001 for all three axes. 

In Fig. 2, the true and estimated yaw angle time histories, and their corresponding estimation 
errors, are shown for a typical run, with and without vector measurement aiding. (The estimated yaw 
angle was computed using the estimated attitude matrix, assuming a 3-2-1 Euler angle sequence). 
Using only GPS measurements, the average yaw estimation error was 7.15 x 10“ 3 deg, with a standard 
deviation of 0.095 deg. When vector measurements were used in combination with the GPS signals, 
the average estimation error was 9.87 x 10” 4 deg, and the estimation error standard deviation 
reduced to 0.022 deg. In Fig. 3, the third component of the angular velocity vector, its estimates 
and corresponding estimation errors are shown for the same run. Using GPS only measurements, 
the steady state estimation error standard deviation was 0.015 deg/s. When vector measurements 
were used in combination with the GPS signals, the estimation error standard deviation reduced to 
0.0065 deg/s (the average rate estimation errors were on the order of 10~ 4 deg/s in both cases). 


Example II 

In this example, the same parameters were used as in Example I, except for the following. The 
three baselines used were now b x = [0.1, 1.0, 0.l] T , 62 = [0.0, 1.0, 0.0] T , 63 = [0.0, 0.0, 1.0] T . As 
can be observed, the first two baselines are almost colinear. The angular rates of the satellite were 
uj = [0, 236,0] deg/hr. The Singer angular acceleration model parameters were set to r = 10 s, 
Cj m — 10 -5 rad/s 2 , pm = po = 001 for all three axes. As in the first example, vector measurements, 
when available, corresponded to a randomly selected, constant vector. 
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**"•(•) T1m» {•) 

(c) Yaw angle (d) Yaw angle estimation error 

Figure 2. Yaw Angle Estimation: (a) and (b) GPS Only Measurements, (c) and (d) 
With Vector Measurement Aiding. 


In Fig. 4, the true yaw angle time history is shown, along with the estimation error time histories 
corresponding to the cases where only GPS measurements were used and where vector observations 
were used along with the GPS measurements. (The estimated yaw angle was computed using the 
estimated attitude matrix, assuming a 3-2-1 Euler angle sequence). As can be observed from Fig. 4, 
the effect of aiding the GPS measurements with vector observations is very substantial in this ill- 
conditioned case. Using only GPS measurements, the average yaw angle steady-state estimation 
error in this run was 7.72 x 10" 3 deg, with an estimation error standard deviation of 0.087 deg. 
When the GPS measurements were aided by vector measurements, the average Euler angle steady- 
state estimation error reduced to 4.6 x 10” 3 deg, with an estimation error standard deviation of 
0.022 deg. In Fig. 5, the estimation error of the third rate component is shown, with and without 
vector observation aiding. Using GPS only measurements, the steady-state rate estimation error 
standard deviation was 9.34 x 10 * deg/s. When vector measurements were used in combination with 
the GPS signals, the standard deviation reduced to 3.51 x 10”^ deg/s (the average rate estimation 
error was on the order of 10~ 5 deg/s in both cases). 
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Figure 3 . 0*3 Estimation: (a) and (b) GPS Only Measurements, (c) and (d) With Vector 
Measurement Aiding. 


CONCLUSIONS 

A nonlinear sequential estimator has been presented, that uses differential GPS carrier phase mea- 
surements to estimate both the attitude matrix and the angular velocity of a spacecraft. The 
algorithm is based on the IRP third -order minimal parametrization of the attitude matrix, which 
is at the heart of its computational efficiency. Avoiding the use of the typically uncertain (and 
frequently unknown) spacecraft dynamic model, the filter uses a polynomial state space model, in 
which the spacecraft angular acceleration is modeled as an exponentially autocorrelated stochastic 
process. When vector observations are available (e.g., from low accuracy Sun sensors or magnetome- 
ters), the estimator’s structure can be easily modified to exploit this additional information and, 
thereby, significantly enhance the algorithm’s robustness and accuracy. Numerical examples have 
been presented, that demonstrate the performance of the proposed algorithm and the advantages of 
aiding the GPS carrier phase signals with vector observations, even when the vector measurements 
are of relatively low accuracy. 
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(a) Yaw angle 




(b) Yaw angle estimation error 


(c) Yaw angle estimation error 


Figure 4. Yaw Angle Estimation: (a) True Angle, (b) GPS Only Measurements, (c) 
With Vector Measurement Aiding. 
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